knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
source('https://raw.githubusercontent.com/oharac/src/master/R/common.R')Read in taxonomic traits filled in by taxon experts and cleaned in prior scripts. Combine with coded sensitivity, adaptive capacity, and exposure from stressor-trait sheets to calculate vulnerability.
_raw_data/xlsx/master_all_taxa_trait_data.xlsx is the raw workbook prepared by Nathalie Butt from the various submissions of the taxa-group experts. This has been processed and cleaned to _data/spp_traits_valid.csv. See earlier scripts in the process.
trait_stressor_rankings/final_scores_all_stressors_traits.xlsx is a workbook with each sheet indicating sensitivity or adaptive capacity; columns in each sheet indicate stressors, and rows indicate traits.
Set up a function to consistently clean trait values. Trait values in the species trait file are already cleaned and adjusted in many cases to get around mismatches; they are generally lower case, no punctuation except for greater/less than signs.
This function also cleans up category and trait names for consistency. All lower case, punctuation and spaces replaced with underscores. The species trait file is already cleaned in this manner.
clean_traitnames <- function(df, overwrite_clean_col = FALSE) {
df <- df %>%
mutate(category = str_replace_all(category, '[^A-Za-z0-9]+', '_') %>% tolower(),
category = str_replace_all(category, '^_|_$', ''),
trait = str_replace_all(trait, '[^A-Za-z0-9]+', '_') %>% tolower(),
trait = str_replace_all(trait, '^_|_$', ''))
if(!overwrite_clean_col & ('trait_value' %in% names(df))) {
return(df) ### without overwriting existing trait_value
}
if(overwrite_clean_col & ('trait_value' %in% names(df))) {
x <- readline(prompt = 'Overwriting existing trait_value column? y/n ')
if(str_detect(x, '^n')) stop('dammit!')
}
### overwrite existing, or add new
df <- df %>%
mutate(trait_value = str_replace_all(tolower(trait_value), '[^0-9a-z<>]', ''))
return(df)
}
clean_traitvals <- function(df) {
x <- df$trait_value
### First: remove numeric commas
y <- str_replace_all(x, '(?<=[0-9]),(?=[0-9])', '') %>%
### then: drop all non-alphanumeric and a few key punctuation:
str_replace_all('[^0-9a-zA-Z<>,;\\-\\.\\(\\)/ ]', '') %>%
### lower case; do it after dropping any weird non-ascii characters:
tolower() %>%
str_trim() %>%
str_replace_all('n/a', 'na') %>%
### convert remaining commas and slashes to semicolons:
str_replace_all('[,/]', ';') %>%
### drop spaces after numbers e.g. 3 mm -> 3mm:
str_replace_all('(?<=[0-9]) ', '') %>%
### drop spaces before or after punctuation (non-alphanumeric):
str_replace_all(' (?=[^a-z0-9\\(])|(?<=[^a-z0-9\\)]) ', '') %>%
### manually fix some valid slashes:
str_replace_all('nearly sessile;sedentary', 'nearly sessile/sedentary') %>%
str_replace_all('live birth;egg care', 'live birth/egg care') %>%
str_replace_all('chitin;caco3mix', 'chitin/caco3 mix') %>%
str_replace_all('0.5-49mm', '0.5mm-49mm')
df$trait_value <- y
return(df)
}
assign_rank_scores <- function(x) {
y <- tolower(as.character(x))
z <- case_when(!is.na(as.numeric(x)) ~ as.numeric(x),
str_detect(y, '^na') ~ NA_real_,
str_detect(y, '^n') ~ 0.00, ### none, NA, no
str_detect(y, '^lo') ~ 0.33,
str_detect(y, '^med') ~ 0.67,
str_detect(y, '^hi') ~ 1.00,
str_detect(y, '^y') ~ 1.00, ### yes
TRUE ~ NA_real_) ### basically NA
return(z)
}Since the species trait file is already cleaned, DO NOT use the clean_traitvals function - it will overwrite the trait_value column.
str_trait_f <- here('trait_stressor_rankings',
'final_scores_all_stressors_traits.xlsx')
str_trait_shts <- readxl::excel_sheets(str_trait_f)sens_traits_raw <- readxl::read_excel(str_trait_f, sheet = 'sensitivity')
sens_traits_df <- sens_traits_raw %>%
janitor::clean_names() %>%
gather(stressor, sens_score, -category, -trait, -trait_value) %>%
mutate(sens_score_orig = as.character(sens_score),
sens_score = assign_rank_scores(sens_score)) %>%
clean_traitnames() %>%
clean_traitvals() %>%
### drop these - going into other categories - we'll change the sheets later
filter(!is.na(category))
str_sens_trait_scores <- sens_traits_df %>%
select(category, trait, trait_value, sens_score, stressor) %>%
mutate(sens_score = ifelse(is.na(sens_score), 0, sens_score)) %>%
filter(!is.na(trait))
### To score for a species/stressor combo, first resolve multiple mutually
### exclusive trait values (using trait_prob) then sum across all traits
spp_sens <- str_sens_trait_scores %>%
left_join(spp_traits, by = c('category', 'trait', 'trait_value')) %>%
filter(!is.na(stressor) & !is.na(taxon)) %>%
group_by(spp_gp, stressor, taxon, trait) %>%
summarize(sens_score = sum(sens_score * trait_prob, na.rm = TRUE)) %>%
group_by(spp_gp, stressor, taxon) %>%
summarize(sens_score = sum(sens_score, na.rm = TRUE)) %>%
ungroup()Unmatched traits between sensitivity scoring sheets and species trait sheets:
x <- spp_traits %>% select(category, trait, trait_value) %>% distinct()
y <- str_sens_trait_scores %>% select(category, trait, trait_value) %>% distinct()These traits are in the species-trait scoring sheets but not found in the sensitivity trait scores (should be adaptive capacity/exposure traits only):
navigation_requirements_sound_or_light_or_magnetic, number_of_sites, number_of_sites_incl_terrestrial_wetlands, adult_mobility, planktonic_larval_duration_pld_exposure, age_to_1st_reproduction_generation_time, are_there_sub_populations, can_the_sex_ratio_be_altered_by_a_stressor, fecundity, global_population_size, lifetime_reproductive_opportunities, max_age, parental_investment, post_birth_hatching_parental_dependence, reproductive_strategy, depth_min_max, eoo_range, zone, if_one_few_size, sub_population_dependence_on_particular_sites
These traits are in the trait-sensitivity scoring sheet but not found in the species scoring (need to be scored for species):
navigation_requirements_sound, navigation_requirements_light, navigation_requirements_magnetic
p <- ggplot(spp_sens, aes(x = stressor, y = sens_score)) +
geom_jitter(size = 1, alpha = .6, height = .1) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5, size = 6)) +
facet_wrap(~ taxon)
ggsave(here('figs/spp_sens_scores.png'), height = 6, width = 6, dpi = 300)
knitr::include_graphics(here('figs/spp_sens_scores.png'))General adaptive capacity traits are basically related to the overall population’s resilience in the face of a threat. Large extents of occurrence, large population sizes, presence of multiple subpopulations, and reproductive strategies fall into this category.
adcap_gen_traits_raw <- readxl::read_excel(str_trait_f, sheet = 'gen_adcap')
adcap_gen_traits <- adcap_gen_traits_raw %>%
select(category, trait, trait_value, adcap_score) %>%
# filter(trait != 'max age' & trait != 'if one/few, size') %>%
mutate(adcap_score_orig = as.character(adcap_score),
adcap_score = assign_rank_scores(adcap_score)) %>%
clean_traitnames() %>%
clean_traitvals()
spp_adcap_gen <- spp_traits %>%
left_join(adcap_gen_traits, by = c('category', 'trait', 'trait_value')) %>%
group_by(spp_gp, taxon) %>%
summarize(adcap_gen_score = sum(adcap_score, na.rm = TRUE)) %>%
ungroup()
# hist(spp_adcap_gen$adcap_gen_score)Unmatched traits between general adaptive capacity scoring sheet and species trait sheets:
x <- spp_traits %>% select(category, trait, trait_value) %>% distinct()
y <- adcap_gen_traits %>% select(category, trait, trait_value) %>% distinct()Traits in species-trait sheets not in general adcap scores: adult_body_mass_body_size, biomineral, calcium_carbonate_structure_location, calcium_carbonate_structure_stages, communication_requirement_sound, extreme_pressure_wave_sensitive_structures, flight, navigation_requirements_sound_or_light_or_magnetic, respiration_structures, adult_mobility, planktonic_larval_duration_pld_exposure, dissolved_oxygen, ph, salinity, sensitivity_to_wave_energy_physical_forcing, thermal_sensitivity_to_heat_spikes_heat_waves, thermal_sensitivity_to_ocean_warming_max_temps_tolerated, can_the_sex_ratio_be_altered_by_a_stressor, feeding_larva_post_hatching_metamorphosis, depth_min_max, zone, air_sea_interface, dependent_interspecific_interactions, extreme_diet_specialization, light_dependence, terrestrial_and_marine_life_stages, if_one_few_size, within_stage_dependent_habitats_condition, across_stage_dependent_habitats_condition
Traits in general adcap scores but not in spp traits:
p <- ggplot(spp_adcap_gen, aes(x = taxon, y = adcap_gen_score)) +
geom_jitter(size = 1, alpha = .6, width = .2, height = .2) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5))
ggsave(here('figs/spp_adcap_gen_scores.png'), height = 6, width = 4, dpi = 300)
knitr::include_graphics(here('figs/spp_adcap_gen_scores.png'))Specific adaptive capacity traits are basically related to an organism’s ability to avoid or mitigate exposure, primarily through movement and larval dispersal.
adcap_spec_traits_raw <- readxl::read_excel(str_trait_f, sheet = 'spec_adcap') %>%
filter(!str_detect(category, 'spatial_scale')) ### drop exposure traits
adcap_spec_traits <- adcap_spec_traits_raw %>%
janitor::clean_names() %>%
gather(stressor, adcap_score, -category, -trait, -trait_value) %>%
mutate(adcap_score_orig = as.character(adcap_score),
adcap_score = assign_rank_scores(adcap_score)) %>%
clean_traitnames() %>%
clean_traitvals() %>%
filter(!is.na(trait))
spp_adcap_spec <- spp_traits %>%
left_join(adcap_spec_traits, by = c('category', 'trait', 'trait_value')) %>%
filter(!is.na(stressor)) %>%
group_by(spp_gp, stressor, taxon) %>%
summarize(adcap_spec_score = sum(adcap_score, na.rm = TRUE)) %>%
ungroup()
# hist(spp_adcap_spec$adcap_spec_score)
# median(spp_adcap_spec$adcap_spec_score); mean(spp_adcap_spec$adcap_spec_score); sd(spp_adcap_spec$adcap_spec_score)
adcap_spec_sum <- spp_adcap_spec %>%
group_by(stressor) %>%
summarize(median = median(adcap_spec_score, na.rm = TRUE),
mean = mean(adcap_spec_score, na.rm = TRUE),
sd = sd(adcap_spec_score, na.rm = TRUE))Unmatched traits between specific adaptive capacity scoring sheet and species trait sheets:
x <- spp_traits %>% select(category, trait, trait_value) %>% distinct()
y <- adcap_spec_traits %>% select(category, trait, trait_value) %>% distinct()Traits in species-trait sheets, not in specific adaptive capacity scores:
adult_body_mass_body_size, biomineral, calcium_carbonate_structure_location, calcium_carbonate_structure_stages, communication_requirement_sound, extreme_pressure_wave_sensitive_structures, flight, navigation_requirements_sound_or_light_or_magnetic, respiration_structures, number_of_sites, number_of_sites_incl_terrestrial_wetlands, dissolved_oxygen, ph, salinity, sensitivity_to_wave_energy_physical_forcing, thermal_sensitivity_to_heat_spikes_heat_waves, thermal_sensitivity_to_ocean_warming_max_temps_tolerated, age_to_1st_reproduction_generation_time, are_there_sub_populations, fecundity, feeding_larva_post_hatching_metamorphosis, global_population_size, lifetime_reproductive_opportunities, max_age, parental_investment, post_birth_hatching_parental_dependence, reproductive_strategy, eoo_range, air_sea_interface, dependent_interspecific_interactions, extreme_diet_specialization, light_dependence, terrestrial_and_marine_life_stages, if_one_few_size, sub_population_dependence_on_particular_sites, within_stage_dependent_habitats_condition, across_stage_dependent_habitats_condition
Traits in specific ad cap scores, not in spp-traits:
p <- ggplot(spp_adcap_spec, aes(x = stressor, y = adcap_spec_score)) +
geom_jitter(size = 1, alpha = .6, width = .2, height = .1) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5, size = 6)) +
facet_wrap( ~ taxon)
ggsave(here('figs/spp_adcap_spec_scores.png'), height = 6, width = 6, dpi = 300)
knitr::include_graphics(here('figs/spp_adcap_spec_scores.png'))| stressor | median | mean | sd |
|---|---|---|---|
| air_temp | 2.34 | 2.079917 | 0.9250551 |
| biomass_removal | 2.00 | 2.325803 | 1.2084672 |
| disease_pathogens | 2.00 | 2.341260 | 1.2125624 |
| entanglement | 2.67 | 2.866005 | 1.1346826 |
| eutrophication_nutrient_pollution | 2.67 | 2.773817 | 1.0897214 |
| habitat_loss_degradation | 3.00 | 3.111772 | 1.3208599 |
| inorganic_pollution | 3.00 | 3.191510 | 1.1682359 |
| invasive_species | 2.00 | 2.129608 | 1.0937445 |
| light_pollution | 3.00 | 2.853341 | 1.1620185 |
| noise_pollution | 3.00 | 3.048002 | 1.3117234 |
| oa | 3.67 | 3.878383 | 1.4350218 |
| oceanographic | 2.67 | 3.032545 | 1.3099011 |
| organic_pollution | 3.00 | 3.191510 | 1.1682359 |
| plastic_pollution | 3.00 | 2.548573 | 1.2221434 |
| poisons_toxins | 3.00 | 3.231926 | 1.1410273 |
| salinity | 3.34 | 3.532366 | 1.2415069 |
| sedimentation | 2.67 | 2.608466 | 1.0894387 |
| slr | 1.00 | 1.449465 | 0.8018596 |
| storm_disturbance | 2.67 | 2.940630 | 1.1640739 |
| uv | 1.67 | 1.613591 | 1.1485770 |
| water_temp | 3.01 | 3.203746 | 1.3396884 |
| wildlife_strike | 1.67 | 1.913639 | 1.1787022 |
Exposure potential modifier checks whether the depth and oceanic zones of the stressor match with the depth and oceanic zones of the species. These fall into the “spatial scale” category with the exception of EOO.
exp_traits_raw <- readxl::read_excel(str_trait_f, sheet = 'spec_adcap') %>%
filter(str_detect(tolower(category), 'spatial.scale')) ### include only exposure traits
exp_traits <- exp_traits_raw %>%
janitor::clean_names() %>%
gather(stressor, exp_score, -category, -trait, -trait_value) %>%
mutate(exp_score_orig = as.character(exp_score),
exp_score = assign_rank_scores(exp_score)) %>%
clean_traitnames()
spp_exposure <- spp_traits %>%
left_join(exp_traits, by = c('category', 'trait', 'trait_value')) %>%
filter(!is.na(stressor)) %>%
group_by(spp_gp, stressor, taxon) %>%
summarize(exposure_mod = as.integer(sum(exp_score, na.rm = TRUE) > 0)) %>%
ungroup() %>%
arrange(stressor, taxon)
non_exposures <- spp_exposure %>%
group_by(taxon, stressor) %>%
mutate(n_gps = n_distinct(spp_gp)) %>%
filter(exposure_mod == 0) %>%
summarize(n_gps_no_exp = n_distinct(spp_gp),
n_gps = first(n_gps),
pct_no_exp = round(n_gps_no_exp / n_gps, 3)) %>%
arrange(stressor, desc(n_gps_no_exp)) %>%
ungroup()Note: this is exposure potential only, based on overlap between species presence and stressor presence - nothing about sensitivity or actual exposure. Check that these logic out.
We will try a calculation for vulnerability \(V\) of species \(i\) to stressor \(j\) that basically looks like this:
\[\text{sensitivity score } S_{i,j} = \mathbf{s}_j^T \mathbf{t}_i\] based on a vector \(\mathbf{s}_j\) of trait-based sensitivity to stressor \(j\), and vector \(\mathbf{t}_i\) of traits of species \(i\);
\[\text{specific adaptive capacity score } K_{i,j} = \mathbf{k}_j^T \mathbf{t}_i\] based on vector \(\mathbf{k}_j\) of trait-based specific adaptive capacity to stressor \(j\); \[\text{general adaptive capacity score } G_{i} = \mathbf{g}^T \mathbf{t}_i\] based on vector \(\mathbf{g}\) of trait-based general adaptive capacity;
\[\text{exposure potential modifier } E_{i,j} = \begin{cases} 1 \text{ when }\mathbf{e}_j^T \mathbf{t}_i > 0\\ 0 \text{ else} \end{cases}\] based on vector \(\mathbf{e}_j\) of trait-based presence of stressor \(j\) (i.e. depth zones and ocean zones in which stressor occurs).
\[\text{vulnerability } V_{i,j} = \frac{S_{i,j} / {S_j}'}{1 + G_i/ {G}' + K_{i,j}/ {K_j}'} \times E_{i,j}\] Each component (\(S_{i,j}, G_i, K_{i,j}\)) is normalized by a reference value (\(S_{j}', G', K_{j}'\) using mean, median, max, etc) for that component for that stressor across all species. Note: median risks referencing to zero for some stressors with few sensitivities (e.g. light pollution); mean risks having a very low reference for the same. Max risks being driven by an outlier, but here the sensitivity scores are generally capped at some low-ish value since there are a finite number of traits that can confer sensitivity. Therefore, we will use max as the reference point. We may wish to consider max possible, which may differ from max observed, in a future iteration?
For species groups with NA in specific adaptive capacity, force to zero (no matching adaptive traits); for species with NA in exposure potential, force to 1 (assume exposure potential).
These results will be saved by species group for now, for future matching to the species level.
### Check that all stressors are matched to ensure proper combining of scores
exp_strs <- spp_exposure$stressor %>% unique()
sens_strs <- spp_sens$stressor %>% unique()
adcap_strs <- spp_adcap_spec$stressor %>% unique()
if(!all(exp_strs %in% sens_strs) | !all(sens_strs %in% exp_strs)) {
stop('Mismatch between stressors in exposure traits and sensitivity traits!')
}
if(!all(adcap_strs %in% sens_strs) | !all(sens_strs %in% adcap_strs)) {
stop('Mismatch between stressors in ad cap traits and sensitivity traits!')
}
if(!all(exp_strs %in% adcap_strs) | !all(adcap_strs %in% exp_strs)) {
stop('Mismatch between stressors in exposure traits and ad cap traits!')
}spp_vulnerability <- spp_sens %>%
left_join(spp_adcap_gen, by = c('taxon', 'spp_gp')) %>%
left_join(spp_adcap_spec, by = c('taxon', 'spp_gp', 'stressor')) %>%
left_join(spp_exposure, by = c('taxon', 'spp_gp', 'stressor')) %>%
### fix NAs
mutate(adcap_spec_score = ifelse(is.na(adcap_spec_score), 0, adcap_spec_score),
exposure_mod = ifelse(is.na(exposure_mod), 1, exposure_mod)) %>%
### calculate means. General adcap is overall; specific adcap and sensitivity are by stressor
group_by(stressor) %>%
mutate(max_adcap_spec = max(adcap_spec_score),
max_sens = max(sens_score)) %>%
ungroup() %>%
mutate(max_adcap_gen = max(adcap_gen_score)) %>%
### rescale components
mutate(sens_rescale = sens_score / max_sens,
adcap_gen_rescale = adcap_gen_score / max_adcap_gen,
adcap_spec_rescale = adcap_spec_score / max_adcap_spec) %>%
### note any stressors with zero for max (e.g. adcap for entanglement),
### will result in a NaN - convert the rescaled score to zero
mutate(sens_rescale = ifelse(is.na(sens_rescale), 0, sens_rescale),
adcap_gen_rescale = ifelse(is.na(adcap_gen_rescale), 0, adcap_gen_rescale),
adcap_spec_rescale = ifelse(is.na(adcap_spec_rescale), 0, adcap_spec_rescale)) %>%
### mash 'em all together
mutate(vuln = sens_rescale / (1 + adcap_gen_rescale + adcap_spec_rescale) * exposure_mod)
x <- spp_vulnerability %>%
filter(taxon == 'reef_fish')
write_csv(spp_vulnerability, here('_output/spp_gp_vulnerability.csv'))vuln_plot <- ggplot(spp_vulnerability, aes(x = stressor, y = vuln)) +
ggtheme_plot() +
geom_jitter(size = 1, alpha = .6, width = .2, height = .02) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5),
strip.background = element_rect(fill = 'grey90')) +
facet_wrap(~ taxon) +
labs(title = 'rescale using max reference point')
ggsave(plot = vuln_plot, filename = here('figs/vuln_plot.png'), width = 8, height = 8, dpi = 300)
knitr::include_graphics(here('figs/vuln_plot.png'))top_3_vuln <- spp_vulnerability %>%
group_by(taxon, stressor) %>%
summarize(vuln = mean(vuln)) %>%
group_by(taxon) %>%
arrange(taxon, desc(vuln)) %>%
filter(vuln >= nth(vuln, 3))
knitr::kable(top_3_vuln)| taxon | stressor | vuln |
|---|---|---|
| cephalopods | eutrophication_nutrient_pollution | 0.3572519 |
| cephalopods | oceanographic | 0.3514481 |
| cephalopods | oa | 0.3340814 |
| corals | light_pollution | 0.5485537 |
| corals | eutrophication_nutrient_pollution | 0.4778458 |
| corals | salinity | 0.4583673 |
| crustacea_arthropods | biomass_removal | 0.3237831 |
| crustacea_arthropods | organic_pollution | 0.2637598 |
| crustacea_arthropods | inorganic_pollution | 0.2619874 |
| echinoderms | water_temp | 0.4404271 |
| echinoderms | eutrophication_nutrient_pollution | 0.4341521 |
| echinoderms | oceanographic | 0.4099282 |
| elasmobranchs | biomass_removal | 0.3548231 |
| elasmobranchs | oceanographic | 0.3204314 |
| elasmobranchs | poisons_toxins | 0.2924205 |
| fish | salinity | 0.3327372 |
| fish | inorganic_pollution | 0.3303307 |
| fish | oa | 0.2978080 |
| marine_mammals | biomass_removal | 0.6092336 |
| marine_mammals | entanglement | 0.5813039 |
| marine_mammals | wildlife_strike | 0.4312690 |
| molluscs | inorganic_pollution | 0.3491484 |
| molluscs | oa | 0.3434182 |
| molluscs | organic_pollution | 0.3419586 |
| plants_algae | light_pollution | 0.4853934 |
| plants_algae | biomass_removal | 0.3094275 |
| plants_algae | organic_pollution | 0.3020574 |
| polychaetes | uv | 0.3969201 |
| polychaetes | habitat_loss_degradation | 0.1552287 |
| polychaetes | plastic_pollution | 0.1367801 |
| reptiles | biomass_removal | 0.5134315 |
| reptiles | invasive_species | 0.4834334 |
| reptiles | entanglement | 0.4576065 |
| seabirds | invasive_species | 0.4145493 |
| seabirds | biomass_removal | 0.3843251 |
| seabirds | slr | 0.3653138 |
| sponges | organic_pollution | 0.2243921 |
| sponges | noise_pollution | 0.2162592 |
| sponges | oceanographic | 0.1848933 |
spp_vuln_summary <- spp_vulnerability %>%
group_by(taxon, stressor) %>%
summarize(n_spp = n_distinct(spp_gp),
sens_mean = mean(sens_rescale),
sens_sd = sd(sens_rescale),
gen_mean = mean(adcap_gen_rescale),
gen_sd = sd(adcap_gen_rescale),
spec_mean = mean(adcap_spec_rescale),
spec_sd = sd(adcap_spec_rescale)) %>%
ungroup() %>%
gather(component1, mean, ends_with('mean')) %>%
gather(component2, sd, ends_with('sd')) %>%
filter(str_sub(component1, 1, 4) == str_sub(component2, 1, 4)) %>%
mutate(component = str_replace(component1, '_mean', '')) %>%
select(-component1, -component2)
ggplot(spp_vuln_summary, aes(x = stressor, color = component, group = component)) +
ggtheme_plot() +
geom_point(aes(y = mean), position = position_dodge(width = .8)) +
geom_linerange(aes(ymin = mean - sd, ymax = mean + sd),
position = position_dodge(width = .8)) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = .5),
strip.background = element_rect(fill = 'grey90')) +
facet_wrap(~ taxon, scales = 'free_y')